% run
% This is a generic parameter file and run script that calls various
% functions to run PARSE reconstructions.  
%
% Change the parameters and file name to change the behavior or inputs.

% clear and close if memory, display are issues
% clear all;
% close all;

recon.Ncpus = 1;            % Number of CPUs available for parallel computing, '1' for serial computing

% characterize imaging experiment
pseq.sw = 103492.884864;    % sampling frequency
pseq.fovcm = 12.8;          % field of view in cm?
pseq.dcoff = 19+9*i;       % DC offset in FID
pseq.offr = offr;             % off-resonance frequency -- units?
pseq.startx = 160;          % ignore this many points at start of FID
pseq.endx = 100;            % ignore this many points at end of FID

% set up image model structure
% reso = 70;                  % image model side dimension
% wfi = 1.3;                  % relative scaling of exponential imaginary
% wfr = 1.3;                  % relative scaling of exponential real
wfi = 1.0;                  % relative scaling of exponential imaginary
wfr = 1.0;                  % relative scaling of exponential real
model = model_setup(reso,pseq.fovcm,wfr+1i*wfi);

% data files
% fid_file = '/home/rishi/data/phantom_studies/Oct2007/s_20071015_01/data/pinecone_01'; %freqoffset=-28;% Foutube Oct 15';
% fid_file = '/home/rishi/data/phantom_studies/Dec2007/s_20071204_02/data/pinecone_01'; % coffee 1.9 gmax
% fid_file = '/home/rishi/data/phantom_studies/Dec2007/s_20071204_01/data/pinecone_01'; % four tube 1.9 gmax
fid_file = '../data/s_20071221_01/data/palm_02'; % four tube 1.9 gmax
trajectory_file = 'KTR1107_19';
% calibration_file = '';    % calibration data in trajectory_file

% parameters for reconstruction algorithm
% recon.alphaM = 10000;     % regularization parameter for magnitude
% recon.alphaf = 5000;      % regularization parameter for frequency map
recon.alphaM = 0;           % regularization parameter for magnitude
recon.alphaf = 0;           % regularization parameter for decay/frequency map
% complex alphaf applies real part to decay and imaginary to freq
recon.wfi = model.wfi;      % preconditioning weight for frequency (?)
recon.wfr = model.wfr;      % preconditioning weight for decay (?)
% recon.NLIST = [(1:35)'*3*120];  % data lengths in progressive-length CG (PLCG)
% recon.FLIST = [.3-((1:35)'/35)*.29];  % thresholds in PLCG


Ta=.06999515;               % duration of acquisition in ms 
delt=1.0/pseq.sw;           % defined from sw in procpar
NNN=floor(Ta/delt);         % number of acquisition points
NNN=NNN-(pseq.startx+pseq.endx);
tim=(0:NNN-1)'*delt;

N1=3*110;
recon.NLIST=[[2*N1:2*N1:18*N1].';NNN];
recon.FLIST=linspace(0.2,0.03,size(recon.NLIST,1)).';

% N1=3*111;
% recon.NLIST=[[N1:N1:13*N1].';NNN];
% recon.FLIST=linspace(0.09,0.04,size(recon.NLIST,1)).';
%recon.NLIST = 3*116*[1;2;3;5;9;17;26;35];  % data lengths in progressive-length CG (PLCG)
%recon.FLIST = [.12;.09;.08;.07;.06;.04;.02;.001];%pretty good.
%recon.FLIST = [.1;.07;.06;.05;.04;.03;.02;.01];
recon.NLENGTHS = size(recon.NLIST,1);  % # of outer loops in PLCG


%recon.FLIST = .02*ones(recon.NLENGTHS,1);
recon.NIT = NIT;            % max # of iterations
recon.GSS = 0;              % use golden section for line search
recon.del = 4e-09;          % step size at which to measure gradient in Cost function CC
recon.bounddecay = 4000/pseq.sw;  % threshold to bound decay parameter
recon.magdown = 0.05;       % reduce mag by this factor if decay exceeded
recon.plot = 0;             % show plots on the fly every recon.plot iterations (0=off)
recon.vis = 2;              % show images on the fly every recon.vis iterations (0=off)
recon.val = 0;              % show criterion values on the fly
recon.time = 4;             % show times on the fly
tic
recon.mat = 'matrix';      % matrix, separable, or filter; save computation or memory
% recon.mat = 'separable';    % matrix, separable, or filter; save computation or memory

% load data and/or calculate phantom FID and trajectory
[pseq.FR,pseq.FI] = LOAD_FID(fid_file);
load(trajectory_file);      % trajectory data and eddy current phases
% phip and phir computed from calibration data using krdpenu19
pseq.krd = real(kss);       % readout k-coordinate
pseq.kpe = imag(kss);       % phase-encode k-coordinate
pseq.phir = phi;            % measured eddy current phase -- r
pseq.phip = 0;              % measured eddy current phase -- p

%% set up useful parameters and preprocess data for reconstruction
frno=1;
fid = fid_setup(pseq,recon.Ncpus,frno);

%% run initialization procedure to get initial estimate if desired
% remember that wf?*param/fid.delt is the decay, frequency, etc.
NN=length(model.ix);
model.pvec = zeros(1,2*NN);
kbasis = make_kbasis(fid,model,recon.mat,recon.Ncpus);
modelsav = model;
%% plcg search
model = cgparse0(fid,recon,model,kbasis);

model = vec2im_model(model);
R2e = real(model.exp)*(-model.wfr/fid.delt);
frmap = imag(model.exp)*(model.wfi/(2*pi*fid.delt));
M0 = abs(model.amp);
M0MSK=(M0>.3*max(max(M0)));

% figure(19); 
% title(strcat('Reconstruction frame : ',num2str(frno)));
% subplot(1,3,1); imagesc(M0); axis image
% subplot(1,3,2); imagesc(R2e); axis image
% subplot(1,3,3); imagesc(frmap); axis image; colormap gray

load handel
sound(y,Fs)

